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Abstract 

^ We address the long-standing problem of computing the region of attraction 

(ROA) of a target set (typically a neighborhood of an equilibrium point) of a con- 
trolled nonlinear system with polynomial dynamics and semialgebraic state and 
input constraints. We show that the ROA can be computed by solving a convex 
linear programming (LP) problem over the space of measures. In turn, this problem 
can be solved approximately via a classical converging hierarchy of convex finite- 
dimensional linear matrix inequalities (LMIs). Our approach is genuinely primal in 
the sense that convexity of the problem of computing the ROA is an outcome of 
optimizing directly over system trajectories. The dual LP on nonnegative continu- 
ous functions (approximated by polynomial sum-of-squares) allows us to generate a 
hierarchy of semialgebraic outer approximations of the ROA at the price of solving a 
^ l sequence of LMI problems with asymptotically vanishing conservatism. This sharply 

l/"") contrasts with the existing literature which follows an exclusively dual Lyapunov 

approach yielding either nonconvex bilinear matrix inequalities or conservative LMI 
conditions. 

OO 
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Given a nonlinear control system, a state-constraint set and a target set (e.g. a neigh- 
borhood of an attracting orbit or an equilibrium point), the constrained controlled region 
of attraction (ROA) is the set of all initial states that can be steered with an admis- 
sible control to the target set without leaving the state-constraint set. The target set 
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can be required to be reached at a given time or at any time before a given timq^J The 
problem of computing the ROA (and variations thereof) lies at the heart of viability the- 
ory (see, e.g., [4J) and goes by many other names, e.g., the reach-avoid or target-hitting 
problem (see, e.g., [21]); in the language of viability theory the ROA itself is sometimes 
referred to as the capture basin 

We show that, in the case of polynomial dynamics, semialgebraic state-constraint, input- 
constraint and target sets, the computation of the ROA boils down to solving an infinite- 
dimensional linear programming (LP) problem in the cone of nonnegative Borel measures. 
Our approach is genuinely primal in the sense that we optimize over state-space system 
trajectories modeled with occupation measures 19]. 

In turn, this LP can be solved approximately by a classical hierarchy of finite-dimensional 
convex linear matrix inequality (LMI) relaxations. The dual LP on nonnegative continu- 
ous functions and its LMI relaxations on polynomial sum-of-squares provide explicitly an 
asymptotically converging sequence of nested semialgebraic outer approximations of the 
ROA. 

Most of the existing literature on ROA computation follows Zubov's approach [22| [32l fT2] 
and uses a dual Lyapunov certificate; see [29j, the survey [TU] . Section 3.4 in [11], and more 
recently [27J and [6j and the references therein. These approaches either enforce convexity 
with conservative LMI conditions (whose conservatism is difficult if not impossible to 
evaluate systematically) or they rely on nonconvex bilinear matrix inequalities (BMIs), 
with all their inherent numerical difficulties. In contrast, we show in this paper that the 
problem of computing the ROA has actually a convex infinite-dimensional LP formulation, 
and that this LP can be solved with a hierarchy of convex finite-dimensional LMIs with 
asymptotically vanishing conservatism. 

We believe that our approach is closer in spirit to set-oriented approaches [7], level-set 
and Hamilton- Jacobi approaches [T71 [2H [23] or transfer operator approaches [30], even 
though we do not discretize w.r.t. time and/or space. In our approach, we model a 
measure with a finite number of its moments, which can be interpreted as a frequency- 
domain discretization (by analogy with Fourier coefficients which are moments w.r.t. the 
unit circle). 

Another way to evaluate the contribution of our paper is to compare it with the recent 
works [T5j [T3] which deal with polynomial approximations of semialgebraic sets. In these 
references, the sets to be approximated are given a priori (as a polynomial sublevel set, 
or as a feasibility region of a polynomial matrix inequality). In contrast, in the current 
paper the set to be approximated (namely the ROA of a nonlinear dynamical system) 
is not known in advance, and our contribution can be understood as an application and 
extension of the techniques of references |TH1 03] to sets defined implicitly by differential 
equations. 

The benefits of our occupation measure approach are overall the convexity of the problem 
of finding the ROA, and the availability of publicly available software to implement and 
solve the hierarchy of LMI relaxations. 

1 The cases of an arbitrarily long but finite time and of asymptotic convergence are not addressed in 
this paper. 
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Our primary focus in this paper is the computation of the constrained finite-time con- 
trolled region of attraction of a given set. This problem is formally stated in Section [2] and 
solved using occupation measures in Section [4| the occupation measures themselves are 
introduced in Section [3) A dual problem on the space of continuous functions is discussed 
in Section [5j The hierarchy of finite-dimensional LMI relaxations of the infinite dimen- 
sional LP is described in Section [6l An extension to the free final time case is sketched in 
Section [T[ Numerical examples are presented in Section [8| and we conclude in Section |9[ 

2 Problem statement 

Consider the control system 

x(t) = fit, x(t),u(t)), x(t) e x, u{t) eu, te [o, T] (l) 

with a given polynomial vector field / with entries f{ G R[i, x, u], i = 1, . . . , n, given final 
time T > 0, and given compact semialgebraic state and input constraints 

x(t) G X := {x G R n : g Xi (x) > 0, i = 1, 2, . . . , n x }, te[0,T] 

u(t) e U := {u e R m : gu^u) > 0, i = 1, 2, . . . , nu }, t e [0,7] U 

with g Xi G R[x], gu i G M[u]. Given a compact semialgebraic target set 

X T := {xeW 1 : g Ti (x) > 0, i = 1, 2, . . . , n T } C X, 

with g Ti G R[x], let 

X(x ) := {x(-) :x(t) = x + / f(r,x(r),u(r))dr, (3) 

Jo 

G [/, x(t) G X, G X T , Vt G [0,T]} 

denote the set of all absolutely continuous admissible controlled trajectories x(-) starting 
from xq, generated by an admissible control u(-) G L 1 ([0, T];M m ). 

The constrained controlled region of attraction (ROA) is then defined as 

X := {x G X : X(x ) + 0}. (4) 

In words, the ROA is the set of all initial conditions for which there exists an admissible 
controlled trajectory. By construction the set X is bounded and unique. 

In the sequel we propose an infinite-dimensional LP approach to computing ROA X 
and show how this reformulation can be approximated by a sequence of LMI problems 
converging to the solution to the LP. 

3 Occupation measures 

In the paper we use the following notations: 
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• Ia(') is the indicator function of a set A, i.e., a function equal to 1 on A and 
elsewhere; 



• A denotes the Lebesgue measure on X C MJ 1 such that 

X(A) — / Ia(x)cI\(x)— I lA{x)dx — / dx 

Jx Jx J A 

is the standard n-dimensional volume of a set A C X; 

• spt fi denotes the support of a measure /i, that is, the closed set of all points x such 
that fi(A) > for every neighborhood A of x. 

3.1 Liouville's equation 

Given an initial condition x and an admissible trajectory x(- | x ) G X(x ) with its 
corresponding control u(- 1 x ) G ^ 1 ([0, T]; R m ) that we assume to be a measurable function 
of #0, define the occupation measure 

fi(A x B x C\x ) := / /AxBxc(^^(^|^o),w(t|x )) ^ 

Jo 

for all subsets i x 5 x C in the Borel cr-algebra of subsets of [0, T] x X x U. Next, for 
a set K let M(K) denote the Banach space of signed Borel measures supported on K 1 
so that a measure v G M(K) can be interpreted as a function that takes any subset of 
K and returns a number in R. Alternatively, elements of M(K) can be interpreted as 
linear functionals acting on the Banach space of continuous functions C{K) 1 that is, as 
elements of the dual space C(K)'. The action of a measure v G M(K) on a test function 
v G C(K) can be modeled with the duality pairing 



(z/, v) :— I v(z) dv(z). 
Jk 



Define further the linear operator C : C^QO, T] x X) C([0, T]xlx[/)by 

i=i 1 

and its adjoint operator £ : C([0, T] x X x U)' — >■ ^([0, T] x X)' by the adjoint relation 
(Ciy,v) :— (v, Cv) — / Cv(t,x,u) dv(t,x,u) 

J[0,T]xXxU 

for all z/ G M([0, T]xXx[/) = C([0, T] x X x [/)' and u G C^flO, T]xl). This operator 
is sometimes expressed symbolically as 



dv ^d[fiiy) dv 

^ c » = -m-2-^x- = -m- dlvf " 

1=1 
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where the derivatives of measures are understood in the sense of distributions (i.e., via 
their action on suitable test functions), and the change of sign comes from the integration 
by parts formula. 

Given a test function v G C 1 ([0,T] x X), it follows from the above definition of the 
occupation measure /i that 

v(T,x(T)) = v(0,x(0))+ [ v(t,x(t\x ))dt 

Jo 

f T 

— v(Q,x(Q)) + / Cv(t, x(t | Xo), u(t | xq)) dt 

Jo 

— v(Q,x(Q))-\- / Cv(t, x : u) d/i(t, x,u\xo). (5) 

J[0,T]xXxU 

Now consider that the initial state is not a single point but that its distribution in space is 
modeled by an initial measure /i G M(X), and that to each initial state x an admissible 
control function u{- \ x ) G L 1 ([0,T];R m ) is assigned in such a way that x(- \ x ) is 
admissibl^J Then we can define the average occupation measure /i G M([0, T] x X x U) 
by 

n(A x B x C) := / fi(A x B x C \ x ) <i/i (^o), (6) 
J x 

and the final measure (It £ M(X T ) by 

/i T (B) := / I B (x(T\xo))dfio(x ). (7) 
Jx 

It follows by integrating ^ with respect to /i that 

/ v(T, x) dfirix) — / t>(0, x) dfio(x) + / Cv^t^x^u) d/i(£,x,ii) 

Jx T Jx J[0,T]xXxU 

or more concisely 

(Ht,v(T,-)) = fa, v(0 r )) + V.e^TlxI), (8) 

which is a linear equation linking the nonnegative measures /i^, /io and /i. Denoting 5 t the 
Dirac measure at a point t and the product of measures, we can write (/i , v(0, •)) = 
(#o /io, and (/ir, v(T,-)) = <8> //t, u )• Then, Eq. Q can be rewritten equivalently 
using the adjoint £ as 

(£'/i, u> = (5 T 0/i T ,v> - <<y ®/Jo, v> Vu G ^([0,^ x X), 

and since this equation is required to hold for all test functions v, we obtain a linear 
operator equation 

C'ji = 8 T <g> fi T - S <g> /io- (9) 



2 The measure fio can be thought of as the probability distribution of Xq although we do not require 
that its mass be normalized to one. 
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This equation is classical in fluid mechanics and statistical physics, where £ is usually 
written using distributional derivatives of measures as remarked above; then the equation 
is referred to as LiouviUe's partial differential equation (PDE). 

Each family of admissible trajectories starting from a given initial distribution /i £ M(X) 
satisfies LiouviUe's equation Q. The converse may not hold in general although for the 
computation of the ROA the two formulations can be considered equivalent, at least from 
a practical viewpoint. Let us briefly elaborate more on this subtle point now. 



3.2 Relaxed ROA 

The control system x(t) — f(t,x(t),u(t)), u(t) G U, can be viewed as a differential 
inclusion 

x(t) G f(t, x(t), U) := {/(*, x(t),u) : ueU}. (10) 

We show in Lemma [4] in Appendix A that any triplet of measures satisfying LiouviUe's 
equation d9|) corresponds to a family of trajectories of the convexified inclusion 



x(t) G cony f(t,x(t),U) (11) 
starting from the initial distribution /i , where conv denotes the convex hull. Let us 



denote the set of absolutely continuous admissible trajectories of (11) by 



X(x ) := {x{-) : x{t) G conv f(t, x{t), U), x{0) = x , x(t) G X x{T) G X T , Vt G [0,T]}. 

Given a famil}|^] of admissible trajectories of the convexified inclusion starting from an 
initial distribution /i , the occupation and final measures can be defined in a complete 
analogy via Q and 0, but now there are only the time and space arguments in the 
occupation measure, not the control argument. In Appendix A, we state and prove the 



correspondence between the convexified inclusion (11) and the measures satisfying the 
Liouville equation pi). 



Define now the relaxed region of attraction as 

X := {xq G X : X(x ) + 0}. 

Clearly X C X and the inclusion can be strict; see Appendix C for concrete examples. 
However, by the Filippov-Wazewski relaxation Theorem [5], the trajectories of the original 



inclusion (10) are dense (w.r.t. the metric of uniform convergence of absolutely continuous 



functions of time) in the set of trajectories of the convexified inclusion (11). This implies 
that the relaxed region of attraction X corresponds to the region of attraction of the 
original system but with infinitesimally dilated constraint sets X and Xt] see Appendix B 
for more details. Therefore, we argue that there is little difference between the two ROAs 
from a practical point of view. Nevertheless, because of this subtle distinction we make 
the following standing assumption in the remainder of the paper. 



Each such family can be described by a measure on C([0, T];M. n ) which is supported on the absolutely 



continuous solutions to (11 ). Note that there may be more than one trajectory corresponding to a single 



initial condition xq since the inclusion (11) may admit multiple solutions 
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Assumption 1 Control system is such that A(X ) — A(X ). 

In other words, the volume of the classical ROA X is assumed to be equal to the volume 
of the relaxed ROA X . Obviously, this is satisfied if X — X , but otherwise these sets 
may differ by a set of zero Lebesgue measure. Any of the following conditions on control 
system is sufficient for Assumption [T] to hold: 

• x(t) G f(t, x(t), U) with f(t, x, U) convex for all t, x, 

• x(t) = f(t,x(t)) + g{t 1 rr(£))w(t), u(t) G U with U convex, 

• uncontrolled dynamics x(t) = f(t,x(t)), 

as well as all controllability assumptions allowing the application of the constrained 
Filippov-Wazewski Theorem; see, e.g., [8] and the discussion around Assumption I in [9j. 



3.3 ROA via optimization 

The problem of computing ROA X can be reformulated as follows: 



q* — sup A(spt /i ) 

s.t. C'fi = S T fi T - 5 O no , 12 x 
fJ> > 0, Mo > 0, ii T > [ } 

spt h C [0, T] x X x U, spt hq C A, spt ht C X Tl 

where the given data are /, A, A T , U, and the supremum is over a vector of nonnegative 



measures (/x, /i , Mt) G M([0, T] x X x [/) x M(X) x M(A T ). Problem ([12]) is an infinite- 
dimensional optimization problem on the cone of nonnegative measures. 



Lemma 1 The optimal value of problem (12) is equal to the volume of the ROA A 0; that 
is, q* = A(A ). 

Proof By definition of the ROA, for any initial condition xq G X there is an admissible 
trajectory in X{x Q ). Therefore for any initial measure /i with spt /i C A there exist an 



occupation measure fi and a final measure fiT such that the constraints of problem (12) 
are satisfied. Thus, q* > A(A ) = A(A ), where the equality follows from AssumptionjlJ 

Now we show that q* < A(A ) = A(A ). Suppose that a triplet of measures (/io,/x, ^t) 
is feasible in (12) and that A(spt/i \ A ) > 0. From Lemma [I] in Appendix A there is 



a family of admissible trajectories of the inclusion (11) starting from /i generating the 
(t, x) -marginal of the occupation measure fi and the final measure (It- However, this is 
a contradiction since no trajectory starting from spt /io \ X can be admissible. Thus, 
A(spt /i \ Xq) = and so A(spt /i ) < A(A )- Consequently, q* < A(A ) = A(A ). □ 
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4 Primal LP on measures 



The key idea behind the presented approach consists in replacing the optimization over 
the support of the measure Ho by the maximization of its mass under the constraint that 
Ho is dominated by the Lebesgue measure. This leads to the following LP: 



P 



sup 
s.t. 



C'h = S T ® Ht - <8> Ho 
H > 0, A > ho > 0, Ht > 

spt h C [0, T] x X x [/, spt /i C X, spt ht C X 5 



(13) 



where the supremum is over a vector of nonnegative measures (/i, /i , /ir) E ([0, T]xlx 



U) x M(X) x M(X T ). In problem ( 13j ) the constraint X > Ho means that A(A) > /io(^) 
for all sets A C X. Note how the objective functions differ in problems (12) and (13). 

The following theorem is then almost immediate. 



Theorem 1 The optimal value of LP problem (13) is equal to the volume of the ROA 
Xq, that is, p* — X(X ). Moreover, the supremum is attained by the restriction of the 
Lebesgue measure to the ROA X . 



Proof: Since the constraint set of problem (13) is tighter than that of problem (12) 



we have by Lemma [I] that A(spt/i ) < A(X ) for any feasible Ho- From the constraint 
Ho ^ A we get Ho{X) — /^o(spt Ho) ^ A(spt/i ) < A(X ) for any feasible Ho- Therefore 
p* < X(Xq). But by definition of the ROA Xq, the restriction of the Lebesgue measure to 
Xq is feasible in (13), and so p* > A(X ). Consequently p* = A(X ). □ 



Now we reformulate problem (13) to an equivalent form more convenient for dualization 



and subsequent theoretical analysis. To this end, let us define the complementary measure 



(a slack variable) p, G M(X) such that the inequality A > Ho ^ in (13) can be written 
equivalently as the constraints Ho + Ao = A, Ho ^ 0, (lq > 0. Then problem (13) is 
equivalent to 



P 



sup 
s.t. 



C'h = S T 



Ho 



Ht - ^o 
Ho + Ho — A 

H > ° 5 > 0, V>t > 0, Ho > 

spt h C [0, T] x X x U, spt Ho C X 1 spt ht C X Tl spt p, c X. 



(14) 



5 Dual LP on functions 



In this section, we derive a dual formulation of problem (14) (and hence (13)) on the 
space of continuous functions. A certain super- level set of any feasible solution to the 
dual problem yields an outer approximation to the ROA Xq. 
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Consider the LP problem 

d* — inf / w(x) d\(x) 
J x 

s.t. Cv(t, x, u) < 0, V (t, x, u) G [0, T] x X x [/ 

>v(0,a;) + l, VxGX ^ 15 ^ 

w(T,x)>0, ViEl T 

> 0, VxGl, 

where the infimum is over (v,w) G C 1 ([0,T] x X) x C(X). The interpretation of the 
dual is intuitive: the constraint Cv < forces v to decrease along trajectories and hence 
necessarily v(0, x) > on X because of the constraint v(T, x) > on X T . Consequently, 
w(x) > 1 on Xq. This instrumental observation is formalized in the following Lemma. 

Lemma 2 If Cv < on [0, T] x X x U, v(T, •) > on X T and w > u(0, •) + 1 on X, 
then w > 1 on X . 

Proof: By definition of X , given any x G X there exists such that x(t) G X, 
G C/ for all t G [0, T] and z(T) G X T . Therefore, since £v<0on [0,T]xIx[/ and 
v(T, •) > on X T , 

< v(T, x{T)) = u(0, x ) + / Cv(t, x{t),u{t)) dt < v(0, x ) < wj(x ) - 1. □ 
We have the following salient result: 



Theorem 2 There is no duality gap between primal LP problems (13) and (14) on mea 



sures and dual LP problem (15) on functions, in the sense that p* — d* . 



Proof: To streamline the exposition, let 

C := C([0,T] x X x U) x C(X) x C(X T ) x C(X), 
M := M([0,T] x X x [/) x M(X) x M(X T ) x M(X), 

and let /C and /C ; denote the positive cones of C and M. respectively. Note that the cone 
K! of nonnegative measures of A4 can be identified with the topological dual of the cone 
/C of nonnegative continuous functions of C. The cone K! is equipped with the weak-* 
topology; see [201 Chapter 5]. 



Then the LP problem (14) can be rewritten as 



p* = sup (7, c) 

s.t. A' j = (3 (16) 



where the infimum is over the vector 

7 := (fi,fi ,fi T ,fL ), 
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the linear operator A' : K! -> ^{[O^T] x X)' x M(X) is defined by 



A' j := (C'n + S <g> /i - S T <8> , Mo + Ao) 



the right hand side of the equality constraint in (16) is the vector of measures 



P := (0, A) e M([0,T] x X) x M(X), 
the vector function in the objective is 

c:= (0,1,0,0) EC, 

and the objective function itself is 



(7, c) = / d/i = fJ>o{X). 
Jx 



The LP problem (16) can be interpreted as a dual to the LP problem 



d* — inf (/?, (v,w)) 

s.t. «4(v, w) — c £ JC. 



(17) 



where the infimum is over (v,w) G C 1 ([0,T] x X) x C(X), and the linear operator 
.A : C^flO, T}xX)x C(X) -> C is defined by 

ytu := {—J£v, w — v(0,-), v(T, •), it;) 



and satisfies the adjoint relation (A/7, u) = (7, .Aw). The LP problem (17) is exactly the 



LP problem (15) 



To conclude the proof we use an argument similar to that of [HI Section C.4]. From 



Theorem 3.10] there is no duality gap between LPs (16) and (17) if the supremum p* 
is finite and the set P := {(^7, (7, c)) : 7 E JC'} is closed in the weak-* topology of 
K! . The fact that p* is finite follows readily from the constraint /i + Ao = A, (lq > 0, 
and from compactness of X. To prove closedness, we first remark that A is weakly-* 
continuous^ since A(v, w) G C for all (v, w) G C 1 ([0, T] x X) x C(X). Then we consider a 
sequence 7^ G and we want to show that its accumulation point \im] £ ^ 00 (A'~fk, (lk, c)) 
belongs to P. Since the supports of the measures are compact and p* is finite (hence 
(i 0j (It and £l are bounded) and since T < 00 (hence fi is bounded), the sequence 7^ is 
bounded and from the weak-* compactness of the unit ball (Alaoglu's Theorem, see [20], 
Chapter 5]) there is a subsequence 7^ that converges weakly-* to an element 7 G Ai so 
that lim^ 00 (^4 / 7 fc ., (7^, c>) G P. □ 



In the above proof we observed that the equivalent LPs (13) and (14) are formulated in 



the dual of a Banach space equipped with the appropriate topology. It follows that the 



supremum in LPs (13) and (14) is attained (by the restriction of the Lebesgue measure 



to X ), a statement already proved by a different means in Theorem [TJ In contrast, the 

4 The weak-* topology on C 1 ([0,T] x X)' x M(X) is induced by the standard topologies on C 1 and 
C - the topology of uniform convergence of the function and its derivative on C 1 and the topology of 
uniform convergence on C. 
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infimum in problem (15) is not attained in C 1 ([0,T] x X) x C(X), but the ^-component 



of feasible solutions of (15) converges to the discontinuous indicator I Xo as we show next. 

Before we state our convergence results, we recall the following types of convergence of a 
sequence of functions Wk '■ X — >- R to a function w : X — >- R on a compact set Id". 
As k —> oo, the functions converge to w: 

• in L 1 norm if J x \wk — w\ dX — >■ 0, 

• in Lebesgue measure if A({x : |itffc(a;) — > e}) — >■ Ve > 0, 

• almost everywhere if 3B C X 1 X(B) = 0, such that Wk — >■ w pointwise on X\B, 

• almost uniformly ifVe > 0,35 C X, A (5) < e, such that Wk —> w uniformly on 
X\B. 

We also recall that convergence in L 1 norm implies convergence in Lebesgue measure and 
that almost uniform convergence implies convergence almost everywhere (see [3l Theorems 
2.5.2 and 2.5.3]). Therefore we will state our results in terms of the stronger notions of 
L 1 norm and almost uniform convergence. 



Theorem 3 There is a sequence of feasible solutions to problem (15) such that its 
w-component converges from above to Ix in L 1 norm and almost uniformly. 

Proof: By Theorem [TJ the optimal solution to the primal is attained by the restriction 
of the Lebesgue measure to X . Consequently, 

P*= [ Ix dX. (18) 
Jx 

By Theorem [2] there is no duality gap (p* = rf*), and therefore there exists a sequence 



(v k ,w k ) e C^flO.Tl x X) x C(X) feasible in (15) such that 



p* — d* = lim / Wk(x)dX(x). (19) 
k ^°° Jx 

From Lemma [2] we have Wk > 1 on X and since Wk > on X, we have Wk > /x on ^ 



for all k. Thus, subtracting (18) from (19) gives 



lim / (wk{x) — Ix {x)) dX(x) — 0, 
k ^°° Jx 



where the integrand is nonnegative. Hence Wk converges to Ix in L 1 norm. From [31 
Theorems 2.5.2 and 2.5.3] there exists a subsequence converging almost uniformly. □ 
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6 LMI relaxations 



In this section we show how the infinite dimensional LP problem ( 14) can be approximated 
by a hierarchy of LMI problems with the approximation error vanishing as the relaxation 
order tends to infinity. The dual LMI problem then yields a converging sequence of outer 
approximations to the ROA. 

The measures in equation (J8J) (or ^) are fully determined by their values on a family 
of functions whose span is dense in C 1 ([0,T] x X). Hence, since all sets are assumed 
compact, the family of test functions in ^ can be restricted to any polynomial basis (since 
polynomials are dense in the space of continuous functions on compact sets equipped with 
the supremum norm). The basis of our choice is the set of all monomials. 

Let MfcN denote the vector space of real multivariate polynomials of total degree less 
than or equal to k. Each polynomial p(x) G can be expressed in the monomial basis 

as 

p{x) = ^2 p a x a , 



\<k 



where a runs over the multi-indices (vectors of nonnegative integers) such that \a\ = 
on < k, and x a = Yli x ai . A polynomial p(x) is identified with its vector of coefficients 
P '■— {Pa) whose entries are indexed by a. Given a vector of real numbers y :— (y a ) 
indexed by a, we define the linear functional L y : R^x] —> R such that 



L y (p) :=p'y := ^2 



PaVc 



where the prime here denotes transposition. When entries of y are moments of a measure 
v, i.e., 

y a = j x a dv(x) 1 

the linear functional models integration of a polynomial w.r.t. v 1 i.e., 

en a 

When this linear functional acts on the square of a polynomial p of degree k 1 it becomes 
a quadratic form in the polynomial coefficients space, and we denote by M^(y) and call 
the moment matrix of order k the matrix of this quadratic form, which is symmetric and 
linear in y: 

L y {p 2 )=p'M k (y)p. 

Finally, given a polynomial g(x) G M[x] we define the localizing matrix Mk(g, y) by the 
equality 

L y {gp 2 ) =p'M k (g, y)p. 
The matrix M^{g 1 y) is also symmetric and linear in y. 

Let y, yo, yr and yo respectively denote the sequences of moments of measures /i, /io, 



and fio such that the constraints in problem (14) are satisfied. Then it follows that these 
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sequences satisfy an infinite-dimensional linear system of equations corresponding to the 
equality constraints of problem (14) written explicitly as 

/ v(T, x) dfirix) — / v(Q, x) dfi (x) — / Cv (t, x, u) d/i(t, x, u) = 0, 

Jx T Jx J[0,T]xXxU 



w(x) dfio(x) + / w(x) dfio(x) = / w(x) d\(x) 
'x Jx Jx 

for the particular choice of test functions v(t, x) = t a x^ and w(x) = x 13 for all a G N and 

j3 G N n . Let us denote by 

A k (y,y ,y Tl y ) = b k 

the finite-dimensional truncation of this system obtained by considering only the test 
functions of total degree less than or equal to 2k. Let further 



d 



Xi 



d, 



Ui 



degg U{ 



dr. 



and consider the primal linear semidefinite programming problem 
p* k = max (2/0)0 



s.t. A k (y,y ,y T ,y ) = b k 












M k (y) h 0, 


M k _ 


- dxt {9xi,y) h 0, 


i = i, 


2,.. 


. ,n x 


M k ^(t(T-t),y) h0, 


M k . 


-duXaui.y) h 0, 


1 = 1 


2,- 


(20) 


M k {y ) h 0, 


M k . 


-d Xl (9xi,yo) h 0, 


z = l, 


2,-. 


. ,n x 


M k {y T ) h 0, 


M k . 


-d Tl {9Ti,yT) h 0, 


i = 1 


2,- 


. ,n T 


M k {y ) h 0, 


M k . 


-d Xi {gxi,yo) h 0, 


z = l, 


2,-. 


.,n x 



where the notation >z stands for positive semidefinite and the minimum is over sequences 
(?/, 7/0, 2/t, 2/o) truncated to degree 2k. Problem (20) is a semidefinite program (SDP), 
where a linear function is minimized subject to convex linear matrix inequality (LMI) 
constraints, or equivalently a finite-dimensional LP in the cone of positive semidefinite 
matrices. 

For the remainder of the section we make the following standard assumption. 

Assumption 2 One of the polynomials modeling the sets X , U resp. X T , is equal to 
g Xi {x) = R 2 X - \\x\\l, gui{u) = Rfj - \\u\\l resp. g Ti {x) = R% - \\x\\\, with R x , Ru resp. 
Rt sufficiently large constants. 

Assumption [2] is made without loss of generality since A, U and X T are bounded, and 
polynomials modeling ball constraints can be added to the semialgebraic descriptions of 
these sets. 



The dual to the SDP problem (20) is given by 

inf w'l 

s.t. —Cv(t, x, u) = p(t, x, u) + q (t, x, u)t(T — t) 

+ ES Qifa x i u )9x l (x) + E£ r *(^ x i u )9Ui{x) 
w(x) - v(0, x) - 1 = p (x) + Y!i=i Qoi{x)g Xi {x) 
v(T, x) = p T (x) + QTi(x)g Ti (x) 
w(x) = s (x) + YJl=i soi(x)g Xi (x), 



(21) 
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where / is the vector of Lebesgue moments over X indexed in the same basis in which 
the polynomial w(x) with coefficients w is expressed. The minimum is over polynomials 
v(t,x) G l^A;^,^] and w G l^fcH, and polynomial sum-of-squares p(t,x,u), qi(t,x,u), 
i = 0, 1, . . . , n x , n(t, x, u), i = 1, . . . , n U: p (x), p T (x), qoi(x), qn{x), s (x), s 0i (x), i = 
1, . . . , nx of appropriate degrees. The constraints that polynomials are sum-of-squares can 
be written explicitly as LMI constraints (see, e.g., [IE]), and the objective is linear in the 



coefficients of the polynomial w(x); therefore problem (21) can be formulated as an SDP 



Theorem 4 There is no duality gap between primal LMI problem (20) and dual LMI 
problem (21), i.e. p\ — d* k . 



In order to prove Theorem [4] we rewrite primal LMI problem (20) in a vectorized form as 
follows 

p* — min c'y 

s.t. Ay = b (22) 
e + Dy G K, 

where y := [y', y' 0l y' Tl y' Q }' and K is a direct product of cones of positive semidefinite ma- 
trices of appropriate dimensions, here corresponding to the moment matrix and localizing 
matrix constraints. The notation e + Dy G K means that vector e + Dy contains entries 
of positive semidefinite moment and localizing matrices, and by construction matrix D 
has full column rank (since a moment matrix is zero if and only if the corresponding 



moment vector is zero). Dual LMI problem (21) then becomes 



d* 



max b'x — e'z 
s.t. A'x + D'z = c 
z G K, 



(23) 



and we want to prove that p* — d* . The following instrumental result is a minor extension 
of a classical lemma of the alternatives for primal LMI (22) and dual LMI (23). The 
notation int K stands for the interior of K. 



Lemma 3 // matrix D has full column rank, exactly one of these statements is true: 

• there exists x and z G int K such that A'x + D'z = c 

• there exists y^O such that Ay = ; Dy G K and c'y < 0. 

Proof of Lemma |3j A classical lemma of alternatives states that if matrix D has full 
column rank, then either there exists z G int K such that D'z = c or there exists y such 
that Dy G K and c'y < 0, but not both, see e.g. [28l Lemma 2] for a standard proof 
based on the geometric form of the Hahn-Banach separation theorem. Our proof then 
follows from restricting this lemma of alternatives to the null-space of matrix A. More 
explicitly, there exists x and z such that A'x + D'z = c if and only if z is such that 
D'z = c with D = DF, c = F'c for F a full-rank matrix such that AF = 0. Matrix D 
has full column rank since it is the restriction of the full column rank matrix D to the 
null-space of A. □. 
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Proof of Theorem [4} First notice that the feasibility set of LMI problem (22) is 
nonempty and bounded. Indeed, a triplet of zero measures is a trivial feasible point 



for (13) and hence (0,0,0, A) is feasible in (14); consequently a concatenation of trun- 
cated moment sequences corresponding to the quadruplet of measures (0, 0, 0, A) is feasible 



in (22) for each relaxation order k. Boundedness of the even components of each moment 



vector follows from the structure of the localizing matrices corresponding to the functions 
from Assumption [2] and from the fact that the masses (zero-th moments) of the measures 
are bounded because of the constraint /i + Ao — A and because T < 00. Boundedness of 
the whole moment vectors then follows since the even moments appear on the diagonal 
of the positive semidefinite moment matrices. 

To complete the proof, we follow [28| Theorem 4] and show that boundedness of the 



feasibility set of LMI problem (22) implies existence of an interior point for LMI problem 



(23), and then from standard SDP duality it follows readily that p* — d* since D has a 
full column rank; see, e.g., [281 Theorem 5] and references therein. 



Let y denote a point in the feasibility set of LMI problem (22), i.e. a vector satisfying 
Ay = b and e + Dy G K. Suppose that there is no interior point for LMI problem (23), 
i.e. there is no x and z G hit K such that A'x + D'z = c. Then from Lemma [3] there 
exists y ^ such that Ay = 0, Dy G K and c'y < 0. It follows that there exists a ray 



y + ty, t > of feasible points for LMI problem (22), hence implying that the feasibility 
set is not bounded. □ 

Now we prove convergence results analogous to those of Theorem [3] as well as a set- wise 



convergence to the ROA X of certain super- level sets of the polynomial solutions to (21 ) 



Theorem 5 Let Wk G M 2 fc[x] denote the w -component of a solution to the dual LMI 



problem (21) and let Wk{x) — min^ Wi{x). Then Wk converges from above to Ix in L 1 



norm and Wk converges from above to Ix in L 1 norm and almost uniformly. 

Proof: From Lemma [2] and Theorem [3J for every e > there exists a (v,w) G C 1 ([0, T] x 
X) x C(X) feasible in (15) such that w > Ix and J x (w — Ix ) d\ < e. Set 



v(t, x) := v(t, x) — et + (T + l)e, 
w(x) := w(x) + (T + 3)e. 



Since v is feasible in (15), we have Cv = Cv — e, and v(T,x) = v(T,x) + e. Since also 



w(x) — £(0,rr) > 1 + 2e, it follows that (v,w) is strictly feasible in (15) with a margin 
at least e. Since [0, T] x X and X are compact, there existj^] polynomials v and w of a 
sufficiently high degree such that \\v — vU^ < e, \\Cv — Cv\\oo < e and \\w — t£>||oo < e. 



The pair of polynomials (v,w) is therefore strictly feasible in (15) and as a result, under 



Assumption 2l feasible in (21) for a sufficiently large relaxation order k (this follows from 



the classical Positivstellensatz by Putinar; see, e.g., [18] or [24]), and moreover w > w. 
Consequently, j x \w — w \ dX < e\(X), and so J x (w — w) dX < eX(X)(T + 4). Therefore 



(w — Ix ) dX < cK 1 w > Ix ; 

Jx 



5 This follows from an extension of the Stone- Weierstrass theorem that allows for a simultaneous 
uniform approximation of a function and its derivatives by a polynomial on a compact set; see, e.g., |16j . 
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where K := [1 + (T + 4)A(X)] < oo is a constant. This proves the first statement since e 
was arbitrary. 

The second statement immediately follows since given a sequence Wk — > Ix i n L\ norm, 
there exists a subsequence that converges almost uniformly to Ix by [31 Theorems 
2.5.2 and 2.5.3] and clearly Wk{x) < min{wfc.(x) : ki < k}. □ 

The following Corollary follows immediately from Theorem [5} 



Corollary 1 The sequence of infima of L MI problems (21) converges monotonically from 



above to the supremum of the LP problem (15), i.e., d* < d* k+1 < d* k and lim^oo d* k = d* . 



Similarly, the sequence of maxima of LMI problems (20) converges monotonically from 
above to the maximum of the LP problem (13), i.e., p* < p* k+1 < p k and lim^oopj = p* . 



Proof: Monotone convergence of the dual optima d* k follows immediately from Theorem [5] 
and from the fact that the higher the relaxation order k, the looser the constraint set of 



the minimization problem (21). To prove convergence of the primal maxima observe that 



from weak SDP duality we have d* k > p k and from Theorems [5] and [2] it follows that 
d k — y d* = p*. In addition, clearly p\ > p* and p* k+1 < p* k since the higher the relaxation 
order k, the tighter the constraint set of the maximization problem (20). Therefore 
Pk ~ * P* monotonically from above. □ 

Theorem [5] establishes a functional convergence of Wk to Ix and Corollary [T] a convergence 
of the primal and dual optima p\ and d k to the volume of the ROA A(X ) = p* — d*. 
Finally, the following theorem establishes a set-wise convergence of the unit super-level 
sets of Wk to Xq. 

Theorem 6 Let Wk G M2fc[#] denote the w -component of a solution to the dual LMI 

(24) 

A(nr = i^ot\^o) = 0. (25) 



problem (21) and let X 0k := {x G 1™ : Wk(x) > 1}. Then X C X 0k , 

hm A(X ofc \X ) = 

and 



Proof: From Lemma |2j we have Wk > Ix and therefore Wk > Ix ok > Ix - From 
Theorem |5j we have Wk — > Ix in L 1 norm on X. Consequently, 

X(Xq) = / I Xo dX = lim / w k dX > lim / Ix ok dX = lim X(X ok ) 

J X k— >oo J j£ k— >oo J j£ k— >oo 

> lim A(njLi*oi) = A(n^° =1 x ofc ). 

k — toe 

But since X C X ok for all k, we must have 

lim A(X ofc ) = X(X ) and A(n^ =1 X fc) = A(X ). 

k— >oo 

This proves the theorem since X ok D X and n^ 1 X A ; D Xq. □ 



16 



7 Free final time 



In this section we outline a straightforward extension of our approach to the problem of 
reaching the target set X T at any time before T < oo (and not necessarily staying in X T 
afterwards) . 

It turns out that the set of all initial states xq from which it is possible to reach Xt at a 
time t <T can be obtained as the support of an optimal solution /1q to the problem 

sup fio(X) 

s.t. C'/i = S T <g> fir - no , . 

V > 0, A > /i > 0, fir > [ } 

spt /iC [0,T] xlx[/, spt fi C X, spt fi T C [0, T] x X T , 

where the supremum is over a vector of nonnegative measures (/i, /i , Mt) G ^([O,^ 1 ] x 



X xU) x M(X) x M([0,T] x X T ). Note that the only difference to problem Q is in 
the support constraints of the final measure \it- 

The dual to this problem reads as 



inf / w(x)d\(x) 



'x 



s.t. Cv(t, x, u) < 0, V (t, x, u) G [0, T] x X x U 

w(x) > v(0,x) + 1, Mx eX 

v(t,x) > 0, V(t,rc) G [0,T] x X T 

> 0, VxGX. 



(27) 



The only difference to problem (15) is in the third constraint which now requires that 
v(t,x) is nonnegative on X T for all t G [0,T]. 

All results from the previous sections hold with proofs being almost verbatim copies. 



8 Numerical examples 

In this section we present three examples to illustrate our approach: a univariate uncon- 
trolled cubic system, an uncontrolled Van der Pol oscillator and a minimum-time control 
of a double integrator. For numerical implementation, one can either use Gloptipoly 3 [Hj 
to formulate the primal problem on measures and then extract the dual solution provided 
by SeDuMi [25J or formulate directly the dual SOS problem using, e.g., YALMIP [31] or 
SOSTOOLS E61 . 



8.1 Univariate cubic dynamics 

Consider the system given by 

x = x(x — 0.5) (a; + 0.5), 
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Figure 1 : Univariate cubic dynamics polynomial approximations (solid line) to the ROA 
indicator function I Xo — ^[-0.5,0.5] (dashed line) for degrees d G {4, 8, 16, 32}. 

the constraint set X — [—1,1], the final time T — 100 and the target set X T — 
[—0.01,0.01]. The ROA can in this case be determined analytically as X = [—0.5,0.5]. 
Polynomial approximations to the ROA for degrees d G {4, 8, 16, 32} are shown in Fig- 
ure [T| As expected the functional convergence of the polynomials to the discontinuous 
indicator function is rather slow; however, the set- wise convergence of the unit super- level 
set is very fast as shown in Table [Tj Note that the volume error is not monotonically 
decreasing indeed what is guaranteed to decrease is the integral of the approximating 
polynomial, not the volume of its unit super- level set. Numerically, a better behavior is 
expected when using alternative polynomial bases (e.g., Chebyshev polynomials) instead 
of the monomials; see the conclusion for a discussion. 

Table 1: Univariate cubic dynamics - relative error of the outer approximation to the ROA 
Xq = [—0.5, 0.5] as a function of the approximating polynomial degree. 



degree 


4 


8 16 


32 


error 


31.60% 


3.32% 0.96% 


1.56% 
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8.2 Van der Pol oscillator 




Figure 2: Van der Pol oscillator polynomial outer approximations (light gray) to the 
ROA (dark gray) for degrees d G {10, 12, 14, 16}. 

As a second example consider a scaled version of the uncontrolled reversed-time Van der Pol 
oscillator given by 

±i = —1x2, 

± 2 = 0.8xi + 10(^1 - 0.21)rr 2 . 

The system has one stable equilibrium at the origin with a bounded region of attraction 

I Cl:= [-1.2, 1.2] 2 . 

In order to compute an outer approximation to this region we take T = 100 and Xt = 
{x : \\x\\2 < 0.01}. Plots of polynomial super-level set approximations of degree d G 
{10, 12, 14, 16} are shown in Figure |2j We observe a relatively fast convergence of the 
super- level sets to the ROA this is confirmed by the relative volume erroij^] summarized 

6 The relative volume error was computed approximately by Monte Carlo integration. 
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Figure 3: Van der Pol oscillator a polynomial approximation of degree 18 of the ROA 
indicator function I Xo . 

in Table [2j Figure [3] then shows the approximating polynomial itself for degree d — 18. 
Here too, a better convergence is expected if instead of monomials, a more appropriate 
polynomial basis is used. 

Table 2: Van der Pol oscillator - relative error of the outer approximation to the ROA Xq as a 
function of the approximating polynomial degree. 



degree 


10 


12 14 


16 


error 


62.7% 


29.6% 20.8% 


14.2% 



8.3 Double integrator 

In our last example we consider a minimum time control of a double integrator 

X\ — x 2 
x 2 = u. 

The goal is to find an approximation of the set of initial states Xq that can be steered to the 
origin in T = 1. Therefore we set X T = {0} and the constraint set such that A C A, e.g., 
X = [—0.7,0.7] x [—1.2, 1.2]. The solution to this problem can be computed analytically 
as 

A = {x : V(x) < 1}, 

where 

ix 2 + 2<Jx~i + \x\ if xi + \x 2 \x 2 \ > 0, 

V( x ) = \ I — — 

I — x 2 + 2-i j—X\ + otherwise. 

The computation results are depicted in Figure [4| again we observe a relatively fast 
convergence of the super- level set approximations, which is confirmed by the relative 
volume errors in Table [U 
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Figure 4: Double integrator polynomial outer approximations (light gray) to the ROA 
(dark gray) for degrees d G {6, 8, 10, 12}. 

9 Conclusion 

The main contributions of this paper can be summarized as follows: 

• contrary to most of the existing systems control literature, we propose a convex 
formulation for the problem of computing the controlled region of attraction; 

• our approach is constructive in the sense that we rely on standard hierarchies of 
finite-dimensional LMI relaxations whose convergence can be guaranteed theoreti- 
cally and for which public-domain interfaces and solvers are available; 

• we deal with polynomial dynamics and semialgebraic input and state constraints, 
therefore covering a broad class of nonlinear control systems; 

• the outer approximation obtained is relatively simple in the sense that it is given 
by a super-level set of a single polynomial of degree given a priori; 

• additional properties (e.g., convexity) of the approximations can be enforced by 
additional constraints on the polynomial (e.g., Hessian being negative definite). 
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Table 3: Double integrator - relative error of the outer approximation to the ROA Xq as a 
function of the approximating polynomial degree. 



degree 


6 8 10 


12 


error 


78.9% 34.6% 21.1% 


16.6% 



The problem of computing the reachability set, i.e. the set of all states that can be 
reached from a given set of initial conditions under input and state constraints, can be 
addressed with the same techniques. Basically, the initial and final measures exchange 
roles. Computation of maximum (controlled) invariant sets should also be amenable to 
our approach. Furthermore, there is a straightforward extension to piecewise polynomial 
dynamics defined over a semialgebraic partition of the state and input spaces one mea- 
sure is then defined for each region of the partition. Our approach should also allow for 
extensions to discrete-time controlled systems, stochastic systems (either discrete-time 
controlled Markov processes or controlled SDEs) and/or uncertain systems. 

The hierarchy of LMI relaxations described in this paper generates a sequence of nested 
outer approximations of the ROA, but it should also be possible, using a similar approach, 
to compute valid inner approximations. 

Numerical examples indicate that the choice of monomials as a dense basis for the set of 
continuous functions on compact sets, while mathematically appropriate (and notationally 
convenient), is not always satisfactory regarding convergence and quality of the approxi- 
mations. However, this is not peculiar to ROA computation problems a similar behavior 
was already observed when computing the volume (and moments) of semialgebraic sets 
in [15] . To achieve better performance, we recommend the use of alternative polynomial 
bases such as (appropriately scaled tensor products of) Chebyshev polynomials. 



Appendix A 

In this Appendix we state and prove the correspondence between the Liouville PDE on 



measures (|9|) and the convexified differential inclusion (11). Let fi(t,x) denote the (t,x)- 
marginal of the occupation measure fi defined through Q, that is, 

p,(A x B) := fi(A x B x U) VAc [0,T], B C X. 

Lemma 4 Let (/i ,A^^t) be a triplet of measures satisfying the Liouville equation 
such that spt /i C X , spt fi C [0, T] x X x U and spt fir C X T . Then there exists a family 



of absolutely continuous admissible trajectories of (11) starting from /i (i.e., trajectories 
in X(xq)) such that the occupation measure and the terminal measure generated by this 
family of trajectories are equal to p, and (It, respectively. 

Proof: First, disintegrate the occupation measure fi as dfi(t,x,u) = dv(u \ t, x)djl{t : x), 
where dv(u \ t, x) is a stochastic kernel, i.e. a probability measure on U for each (t, x) 
which is measurable in (t, x) for a fixed first argument. Then we can rewrite equation Q8J) 
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as 



v(T,') dfiT 



A x 



v(0, •) dp + 
v(0, •) dp + 
v(0, •) dp + 



+ grad v • /(£, x, u) dv{u \ t, x) dp(t, x) 



'x 



[0,T]xX 



dv 
di 
dv 
di 



+ grad v 



f(t,x,u)dis(u\t,x) dp(t,x) 



+ grad v • f(t, x) dp(t, x) , 



(28) 



where 



f(t,x):— / f(t, x, u) dv{u \t,x) G conv f(t, x, U). 
Ju 

Therefore we will study the trajectories of the differential equation 

x(t) = f(t,x(t)). 



(29) 



In the remainder of the proof we show that the measures \±t and p are generated by a 
family of absolutely continuous trajectories of this differential equation (which is clearly 
a subset of trajectories of the convexified inclusion (11)) starting from p . Note that the 
vector field / is only known to be measurable, so this equation may not admit a unique 
solution. 

Observe that the ^-marginal of p (and hence of p) is equal to the Lebesgue measure 
restricted to [0, T] scaled by p := po(X) (=fir(X)) - this follows by plugging in the family 
of test functions v(t,x) = t k , k G N, into equation Q. Therefore we can disintegrate p 
as 

dp,(t,x) = dp t (x)dt, (30) 

where dp t (x) is a stochastic kernel on X given t scaled by p and dt is the standard 
Lebesgue measure. The kernel p t can be thought of as the distribution^ of the state at 
time t. The kernel p t is defined uniquely dt-almost everywhere, and we will show that 
there is a version such that the function t \-> J x w(x) dp, t (x) is absolutely continuous for 
all w G C X {X) and such that the continuity equation 

— f w(x)dp t (x) = [ gTdidw(x) • f(t, x) dp t (x) Vw G C l (X) (31) 
dtJx Jx 

is satisfied almost everywhere on [0,T] with the initial condition p . 

Fix w G C 1 (X) and define the test function v(t,x) := ijj(t)w(x), where ijj G C 1 ([0,T]). 
Then from equation (28) 



d(%l)w) 



dt 



+ grad(^w) • f(t, x) dp(t, x) 



i/j(T)w dfiT — / i/j(0)w dpo + 

X T JX J[0,T]xX 

= ijj(0) / wdpQ+ I I ip(t)w(x) + ip(t)gTa,dw(x) • f(t, x) dp t {x)dt 
Jx Jo Jx 

T I 

dt, 



■0(0) / w dpo + 
'x Ju 



ip I w dp t + ip grad w • / dp t 
x Jx 



7 It will become clear from the following discussion that for t — and t — T this kernel (or a version 
thereof) coincides with /iq and /x^, respectively; hence there is no ambiguity in notation. 
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which can be seen as an equation of the form 



^(T)d = <0(O)c+ [ Tp(t)a(t)+i>{t)b(t)dt W> e C^QO, T\), (32) 
Jo 

where c := J x w(x)dfio(x), d :— J Xt w(x) dfiT and b(t) :— f x gr&dw • f(t,x)dfi t (x) are 
constants and a(t) is an unknown function. One solution is clearly a(t) = J x wdfi t . Now 
we show that 



a(t):=c + / b(r) dr = / w dfi + / / gr ad w •/ d\i T dr 
Jo Jx Jo Jx 

also solves the equation. Indeed, since from (28) with v replaced by w we have a(T) — 
J x w dfir — d, integration by parts gives 

T rT 

ip(t)a(t) dt = i/j(T)d - i/j(0)c - / i/j(t)b(t) dt, 



so a(t) indeed solves equation (32). Now we prove that this solution is unique. Since a is 
a solution we have 

^(T)d = ^(0)c+ / ip(t)a{t)+ip(t)b(t)dt, 



and subtracting this from (32) we get 



or equivalently 



0= f ij;(t)[a{t)-a{t)]dt V^G^ftT]), 
Jo 

0= f cf)(t)[a(t)-a(t)]dt V0EC([O,T]). 
Jo 



Since C([0, T]) is dense in L 1 ([0, T]), this implies a(t) = a(t) ^-almost everywhere. Con- 
sequently, since is separable, 

/ w dfi t = / w(x)d/iQ + / / gr&dw • f dfi T dr \/wEC 1 {X) (33) 

Jo Jx 

dt-almost everywhere. The right-hand side of this equality is an absolutely continuous 
function of time for each w G C 1 (X) and the left-hand side is a bounded positive linear 
functional on C(X) for all t G [0, T\. By continuity, the right-hand side is a bounded 
positive linear functional on C 1 (X) for all t G [0, T] which can be uniquely extended to 
a bounded positive linear function on C(X) (since C 1 is dense in C). Therefore, for all 
t G [0, T] the right-hand side has a representing measure and hence there is a version of 
fit such that the equality (33) holds for all t G [0, T\. With this version of fi t the function 
t i — ^ J x w(x) dfi t (x) is absolutely continuous and fi t solves the continuity equation (31). 

To finish the proof, we use [H Theorem 3.2] which asserts the existence of a nonnega- 
tive measure a on C([0,T];R n ) which corresponds to a family of absolutely continuous 
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solutions to ODE (29) whose projection at each time t G [0,T] coincides with \i t . More 



precisely, there is a nonnegative measure a G M(C([0, T\\ W 1 )) supported on a family of 



absolutely continuous solutions to ODE (29) such that for all measurable w : MJ 1 —> 



/ w{x)nt{x)= I w{x{t))da{x{-)) V*e[0,T]. (34) 

J X Jc([0,T];R n ) 



It follows from (30) that 



fi(A x B) — I lA(t)lB(x) d/i(t, x) — I / Ib{x) dfi t (x) dt. 

J[0,T]xX JO JX 



Therefore, using (|34|) with w(x) = Ib(x) and Fubini's theorem, we get 

'C([0,T] ; 



fi(AxB)= f f I AxB (t,x(t))dtda(x(-)). 

JC(\0,T]:R n ) JO 



and so the occupation measure of the family of trajectories coincides with Jl. Clearly, 
the initial and the final measures of this family coincide with /i and \±t as well. As a 
result cr-almost all trajectories of this family are admissible. The proof is completed by 
discarding the null-set of trajectories that are not admissible, which does not change the 
measure a and the generated measures p,, /io, ^t- D 



Appendix B 



In this Appendix we elaborate further on the discussion from Section 3.2 on the connection 
between the classical ROA and the relaxed ROA. Let us recall the definition of the classical 
ROA 

X := {xq G X : X(x ) + 0}, 

where 

X(x ) := {x(-) : x{t) G f(t, x(t), U), x(0) = x , x(t) G X, x{T) G X T }, 

and x(-) is required to be absolutely continuous. Similarly, recall the definition of the 
relaxed ROA 

X := {xq G X : X(x ) + 0}, 

where 

X(x ) := {x(-) : x(t) G conv f(t, x(t), U), x(0) = x , x(t) G X, x(T) G X T } 
with x(-) absolutely continuous. Obviously, it holds 

X C X , (35) 
and the question is whether this inclusion is strict or not. 
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Denote B e (a) := {x G MJ 1 : \\x — a\\ 2 < e} and define the dilated constraint sets 

r:=l0 5 e (O) and I^:=I T ©5 £ (0), 

where denotes the Minkowski sum of two sets. Accordingly, the dilated ROA and the 
dilated relaxed ROA are 

X e := {x G X : X'(x Q ) + 0}, 
X e := {x G X : X*(x ) + 0}, 

where 

X e (x ) := {x{-) : x(t) G /(t, x{t), U), x(0) = x , x(t) G X e , x(T) G Xf}, 
^(xq) := {x(-) : G conv / (t, x(t), U), x{0) = x , x{t) G X e , G X^}. 



Since the constraint sets are compact and the vector field / Lipschitz, it follows from the 



equivalence between the trajectories of the convexified inclusion (11) and solutions to the 
Liouville equation Q, stated in Lemma [4] of Appendix A, and from Filippov-Wazewski's 
relaxation Theorem (see, e.g., [5]) that 

Xq = spt/io C P|^o- 

e>0 

In contrast, for all e > it holds 



In general inclusion (35) is strict. However, we argue that for most practical purposes 
the relaxed ROA X and the true ROA X are the same. Indeed, for any x G X there 
exists a sequence of admissible control functions Uk(-) such that 

sup distx{xk{t)) — > and distx T (xk(T)) — >■ as k — >• oo, 
te[o,T] 

where Xk(-) denotes the solution to the ODE ([!]) corresponding to the control function 
Ufc(-), and dist^(^) := inf { || ^ — x\\ 2 : 2; G ^4} denotes the distance to a set A. 



Appendix C 

In this Appendix we describe two contrived examples of control systems ([!]) for which 
the relaxed ROA X is strictly larger than the classical ROA X ; see Appendix B for 
definitions. 

Let f(t,x,u) — u, U — { — 1, +1}, X — X T = {0} for, e.g., T = 1. Obviously there is no 
admissible trajectory in X(0), whereas there is a feasible triplet of measures satisfying Q 
given by /i = ^0, Mr = ^0 and fi = A[ ,i] So + £+1), where A[ ,i] denotes the 

restriction of the Lebesgue measure to [0,1]. Therefore in this case X = ^ X = 
{0}, but A(X ) = A(X ). Assumption [l] is therefore satisfied. Note that the relaxed 
solution corresponds to an infinitely fast chattering of the control input between —1 and 
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+1 which can be arbitrarily closely approximated by chattering solutions of finite speed; 
the singleton constraint set X 1 however, renders such solutions infeasible. 

Another example for which the gap (e.g., in volume) between X and X can be as large 
as desired is the following. Consider x — u GR 2 with 

xeX := B R ((-1 - R, 0)) U [-1, 1] x {0} U #i((+2, 0)) C R 2 

for a given R > and with u G U := {-1, l} 2 and X T := Bi((+2, 0)). Then X = X T is 
strictly smaller than X — X, and A(X ) = ir, whereas X(Xq) — (1 + R 2 )tt. Assumption [l] 
is therefore not satisfied for R > 0. In this example, regular solutions starting in the 
left ball cannot transverse the line to the right ball; this is, by contrast, possible for the 
relaxed solutions using an infinitely fast chattering. 
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